library(magrittr)
library(stringr)
library(plyr)
library(tidyverse)
library(readxl)
library(stargazer)
library(multiwayvcov)
library(lmtest)
library(Amelia)


rm(list=ls())
home = 'C:/Users/Jason/Dropbox/VNA_Responsiveness/Analysis/JOP-dataverse/'


survey = paste0(home, 'survey-outcomes.xlsx') %>%
  read_xlsx %>%
  mutate(Citizen=mapvalues(x=Treatment, 
                           from=c('Control',
                                  'Citizen',
                                  'Firm'), 
                           to=c(0,
                                1,
                                0)),
         Citizen=as.integer(Citizen),
         Firm=mapvalues(x=Treatment, 
                        from=c('Control',
                               'Citizen',
                               'Firm'), 
                        to=c(0,
                             0,
                             1)),
         Firm=as.integer(Firm),
         Treatment=factor(x=Treatment, 
                          levels=c('Control',
                                   'Citizen',
                                   'Firm')))


# Columns 1-3
fit_1 = lm(Q1 ~ Citizen + Firm, 
           data=subset(survey, 
                       !is.na(Treatment)))
fit_2 = lm(Q1 ~ Citizen + Firm + FullTime + CentNom + Competitive, 
           data=subset(survey, 
                       !is.na(Treatment)))
fit_3 = lm(Q1 ~ Citizen*Prop.Citizen + Firm*Prop.Firm + FullTime + CentNom + Competitive, 
           data=subset(survey, 
                       !is.na(Treatment)))
fit_1_cluster = cluster.vcov(model=fit_1, ~Province)
fit_2_cluster = cluster.vcov(model=fit_2, ~Province)
fit_3_cluster = cluster.vcov(model=fit_3, ~Province)


# Column 4: multiple imputation with central nominees
set.seed(31415)
mi_data_4 = subset(survey, 
                   !is.na(Treatment),
                   select=c(ID,
                            Q1,
                            Citizen,
                            Firm,
                            FullTime,
                            CentNom,
                            Competitive,
                            Prop.Citizen,
                            Prop.Firm)) %>%
  amelia(m=100, 
         idvars='ID', 
         noms=c('Q1',
                'Citizen',
                'Firm',
                'FullTime',
                'CentNom',
                'Competitive'), 
         p2s=0)
mi_fits_4 = llply(.data=mi_data_4$imputations, 
                  .fun=function(x) {
                    fit_data = merge(x=x, 
                                     y=subset(survey,
                                              select=c(ID,
                                                       Province)), 
                                     by='ID', 
                                     all.x=T)
                    fit = lm(Q1 ~ Citizen*Prop.Citizen + Firm*Prop.Firm + FullTime + CentNom + Competitive, 
                             data=fit_data)
                    fit_cluster = cluster.vcov(model=fit, ~Province)
                    coeftest(x=fit, vcov.=fit_cluster)
                  }, .progress='text', .inform=T)
mi_melded_4 = mi.meld(q=ldply(.data=mi_fits_4, 
                              .fun=function(x) x[,'Estimate'] )[,-1],
                      se=ldply(.data=mi_fits_4, 
                               .fun=function(x) x[,'Std. Error'] )[,-1]) %>%
  do.call(rbind, .) %>%
  t %>%
  as.data.frame %>%
  set_colnames(c('Beta',
                 'SE')) %>%
  mutate(tValue=Beta/SE,
         pValue=pt(q=abs(tValue), 
                   df=length(coef(object=fit_3))-1, 
                   lower.tail=F)) %>%
  sapply(round, digits=3) %>%
  set_rownames(names(coef(object=fit_3)))


# Columns 5-7
fit_5 = lm(Q1 ~ Citizen + Firm, 
           data=subset(survey, 
                       !is.na(Treatment) & CentNom==0))
fit_6 = lm(Q1 ~ Citizen + Firm + FullTime + CentNom + Competitive, 
           data=subset(survey, 
                       !is.na(Treatment) & CentNom==0))
fit_7 = lm(Q1 ~ Citizen*Prop.Citizen + Firm*Prop.Firm + FullTime + CentNom + Competitive, 
           data=subset(survey, 
                       !is.na(Treatment) & CentNom==0))
fit_5_cluster = cluster.vcov(model=fit_5, ~Province)
fit_6_cluster = cluster.vcov(model=fit_6, ~Province)
fit_7_cluster = cluster.vcov(model=fit_7, ~Province)


# Column 8: multiple imputation withOUT central nominees
set.seed(31415)
mi_data_8 = subset(survey, 
                 !is.na(Treatment) & CentNom==0,
                 select=c(ID,
                          Q1,
                          Citizen,
                          Firm,
                          FullTime,
                          Competitive,
                          Prop.Citizen,
                          Prop.Firm)) %>%
  amelia(m=100, 
         idvars='ID', 
         noms=c('Q1',
                'Citizen',
                'Firm',
                'FullTime',
                'Competitive'), 
         p2s=0)
mi_fits_8 = llply(.data=mi_data_8$imputations, 
                .fun=function(x) {
                  fit_data = merge(x=x, 
                                   y=subset(survey,
                                            select=c(ID,
                                                     Province)), 
                                   by='ID', 
                                   all.x=T)
                  fit = lm(Q1 ~ Citizen*Prop.Citizen + Firm*Prop.Firm + FullTime + Competitive, 
                           data=fit_data)
                  fit_cluster = cluster.vcov(model=fit, ~Province)
                  coeftest(x=fit, vcov.=fit_cluster)
                }, .progress='text', .inform=T)
mi_melded_8 = mi.meld(q=ldply(.data=mi_fits_8, 
                            .fun=function(x) x[,'Estimate'] )[,-1],
                    se=ldply(.data=mi_fits_8, 
                             .fun=function(x) x[,'Std. Error'] )[,-1]) %>%
  do.call(rbind, .) %>%
  t %>%
  as.data.frame %>%
  set_colnames(c('Beta',
                 'SE')) %>%
  mutate(tValue=Beta/SE,
         pValue=pt(q=abs(tValue), 
                   df=length(coef(object=fit_3))-2, 
                   lower.tail=F)) %>%
  sapply(round, digits=3) %>%
  set_rownames(setdiff(names(coef(object=fit_3)),
                       'CentNom'))


# Write LaTeX table -- with Columns 4 and 8 blank -- to clipboard
stargazer(coeftest(x=fit_1, vcov.=fit_1_cluster),
          coeftest(x=fit_2, vcov.=fit_2_cluster),
          coeftest(x=fit_3, vcov.=fit_3_cluster),
          coeftest(x=fit_3, vcov.=fit_3_cluster),
          coeftest(x=fit_5, vcov.=fit_5_cluster),
          coeftest(x=fit_6, vcov.=fit_6_cluster),
          coeftest(x=fit_7, vcov.=fit_7_cluster),
          coeftest(x=fit_7, vcov.=fit_7_cluster),
          title='Debate preparation. Forms the basis of Figures 5, 8, and 9.', 
          covariate.labels=c('Citizen',
                             '% Citizen',
                             'Citizen $times$ % Citizen',
                             'Firm',
                             '% Firm',
                             'Firm $times$ % Firm',
                             'Full-time',
                             'Central Nominee',
                             'Competitive',
                             'Constant'),
          order=c(1:2, 8, 3:4, 9, 5:7, 10),
          coef=list(NULL, NULL, NULL, NA, NULL, NULL, NULL, NA),
          se=list(NULL, NULL, NULL, NA, NULL, NULL, NULL, NA),
          align=T, 
          header=F,
          no.space=T, 
          star.cutoffs=NA,
          column.labels=c('All delegates','No central nominees'),
          column.separate=c(4, 4),
          dep.var.labels.include=F,
          label='appendix-table-05.1',
          notes.append=F, 
          notes='Province-clustered standard errors in parentheses;') %>%
  writeClipboard


# Columns 4 and 8 must be manually entered
mi_melded_4[c('Citizen',
              'Prop.Citizen',
              'Citizen:Prop.Citizen',
              'Firm',
              'Prop.Firm',
              'Firm:Prop.Firm',
              'FullTime',
              'CentNom',
              'Competitive',
              '(Intercept)'),c('Beta','SE')]
rbind(mi_melded_8[c('Citizen',
                    'Prop.Citizen',
                    'Citizen:Prop.Citizen',
                    'Firm',
                    'Prop.Firm',
                    'Firm:Prop.Firm',
                    'FullTime'),c('Beta','SE')],
      c(NA, NA), # placeholder for missing CentNom coefficient
      mi_melded_8[c('Competitive',
                    '(Intercept)'),c('Beta','SE')])
